OPTIMUM  DIGITAL  CONTROL  SYNTHESIS 


•k    k    A  A  A  & 


Captain  Don  J.  Ogden  USMC 


OPTIMUM  DIGITAL  CONTROL  SYNTHESIS 


By 


Don  J.  Ogden 

// 
Captain,  United  States  Marine  Corps 


Submitted  in  partial  fulfillment  of 
the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE 

IN 

ENGINEERING  ELECTRONICS 

United  States  Naval  Postgraduate  School 
Monterey,  California 

19  6  5 


Li  brary 

U.  S.  Naval  Postgraduate  SciuMl 

Monterey,  California 


OPTIMUM  DIGITAL  CONTROL  SYNTHESIS 
by 
Captain  Don  J.  Ogden  USMC 

This  work  is  accepted  as  fulfilling 
the  thesis  requirements  for  the  degree  of 
MASTER  OF  SCIENCE 
IN 
ENGINEERING  ELECTRONICS 
from  the 
United  States  Naval  Postgraduate  School 


ABSTRACT 
Sampled-data  control  systems  are  coming  of  age.   Many  methods 
have  been  formulated  by  which  both  continuous  and  sampled-data 
control  systems  may  be  optimized.   Here  one  method  of  optimum 
control  design  is  presented  using  the  approach  of  dynamic  pro- 
gramming.  A  program  to  generate  the  state  transition  matrices 
from  the  system  dynamics  was  developed.   A  program  for  minimizing 
rather  general  cost  functions  was  written  and  examples  are  pre- 
sented utilizing  the  results. 
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CHAPTER  I 
INTRODUCTION 

In  general,  any  linear  control  system  may  be  expressed  by  the  fol- 
lowing state  equation: 

x  =  Fx  +  Du  (1) 

where  x.  represents  the  system  variables  and  u.  are  the  control  vari- 
ables. F  and  D  are  nxn  and  nxm  matrices  defining  the  dynamical  rela- 
tionship of  the  states  and  controls.  We  will  sample  the  states  at  dis- 
crete intervals  t,  and  at  a  sampling  rate  with  period  T.  A  recursive 
relationship  from  the  solution  of  the  above  differential  equation  may 
be  established: 

x(k+l)=  eFTx(k)  +  fT  eF(T"T)Du(k)dT  (2) 


where  p2  T2         pn  -j.n 

eFT=  I  +  FT  +  —  +•  •  •+  T"  (3> 

2!  n! 


Equation  (2)  may  be  written: 

x(k+l)=  (j>x(k)  +  Au(k)  (4) 


FT 
where  <j>  =  e  (5) 


/T  f(T-t) 
and  A  -  (   e      Ddx  (6) 


As  an  example  of  the  above  equations,  the  open-loop  system  of 
Figure  1-1  will  be  considered. 


i/s  L 


Figure  1-1:  Block  diagram  of  an  open  loop  control  system. 
By  inspection,  the  following  state  equation  defines  the  example  sys- 


tem: 


x  = 


0 

1 

x  + 

0 

0 

0 

_l_ 

(7) 


The   <J>  matrix  as  defined  by  (5)  has  only  two  non-zero  terms  in  the 
Taylor  expansion  for  the  exponential  function: 


cf,  = 


1   0 
0   1 


0   1 
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0 


fil 


(8) 


or 


<j>  = 


1    T 
0   1 

In  a  like  manner,  the  convolution  integral  of  (6)  may  be  solved 
for   A  : 


(9) 


A  = 


giving 


1 

T-T 

0 

p      1  _ 

1 

'2 

T 

dx 


(10) 


(ID 


The  system  difference  equation  may  now  be  written  as  a  function 
of  the  sampling  period  T: 


x(k+l)  = 


1   T 
0   1 


x(k)  + 


u(k) 


(12) 


The  calculation  of  the  transition  matrix,  <}>,  and  the  external 
forcing  function's  impulse  response  matrix,  A,  from  the  matrices  F  and  D, 


of  the  process  dynamics  was  accomplished  with  the  digital  computer.   This 
program  required  an  input  of  F  and  the  sampling  rate  T.   The  amount  of  ac- 
curacy desired  may  be  varied  by  the  user  by  inputing  the  least  significant 
figure  desired.   This  test  point  for  accuracy  results  from  the  fact  that 
the  sum  of  all  terms  in  a  Taylor  series  following  the  test  term  is  less 
in  magnitude  than  the  test  term.    Figure  1-2  is  the  flow  chart  of  the 
program.   The  program  has  been  incorporated  into  other  programs  written 
and  presented  in  this  paper. 


START 

Y 

READ: 

N,DT,TEST,F,  and  D 


TM=  0.0 
TERM=  I 
WORM=  0.0 
<j)  =  TERM 


TELM=  TERMED T 
DELP=  TELM 
DELM=  0.0 


TM=   TM+1.0 

DELM=  DELM  -   TELM*F*DT/ (TM+1 . 0) 

WORM= WORM+T ERM* F *DT / TM 


TERM=WORM 
WORM=  0.0 
TELM=DELM 


DELP=DELP+TELM 
DELM=   0.0 
<t>=<jH-TERM 


ABC   =    |TERM(IR,IC)lmax 


f       TEST: ABC     ) 


A=   A+4>*DELP*D 


PRINT:      <J>   and   A 


A 

STOP 


Figure   1-2:      Flow  Chart   for  PROGRAM  PHIDEL 


The  optimal  control  of  a  system  depends  upon  a  valid  determination 
of  an  optimal  control  policy  which  will  maximize  or  minimize  a  perform- 
ance criteria.   This  control  policy  may  be  generated  from  criteria  in- 
volving time,  system  state  variables,  noise,  control  effort,  etc.   In 
some  cases,  this  control  policy  may  be  generated  by  active  networks  or 
filters  incorporated  within  the  system,  or  a  digital  computer  may  be 
used  to  generate  the  control.   To  accomplish  this,  all  dynamics  of  the 
system  must  be  fed  back  to  the  digital  computer  as  indicated  in  Figure  1-3 


Figure  1-3:  Multivariable  Digital  Control  System 
A  general  approach  to  the  optimization  of  any  dynamic  system  was 
developed  through  the  efforts  of  R.  E.  Bellman  and  has  been  called 
"dynamic  programming".   As  an  example  of  this  method,  in  a  missile  prob- 
lem, one  might  choose  a  proper  cost  function  (terminal  CEP)  and  minimize 
this  cost  throughout  the  trajectory.   The  optimum  trajectory,  according 
to  the  chosen  cost  function,  would  be  computed. 

Dynamic  programming  appears  to  be  an  excellent  method  for  determ- 
ining an  optimal  control  process  because  of  its  adaptability  to  use  on 
a  digital  computer.   Throughout  the  remainder  of  this  paper,  this  tech- 
nique will  be  utilized  in  formulating  general  programs  to  synthesize 
optimum  controllers. 


CHAPTER     II  ' 
PROGRAM  OPCON1 

Again  let  us  consider  the  system: 

x  =  F  x  +  D  u  ( 1 '  ) 

whose  discrete  solution  is 

x(k+l)=  <j>x(k)  +  Au(k)  (4') 

where  we  now  take  the  case  where  u(k)  is  a  scalar  control  and  A_  is 
a  column  vector  of  system  impulse  responses  due  to  control  inputs. 
The  recursion  formulas  will  be  derived  for  this  system  minimizing  the 
sum  of  the  quadratic  terms  as  described  by  the  matrix  Q  over  an  N  stage 
process.   The  Q  matrix  is  specified  by  the  requirements  of  the  user  and 
the  application  of  solid  engineering  judgment.   This  cost  function  may 
be  written: 

N   T 
J(N)=  minimum   I   x i(k)Qx(k)  (13) 

k=l  " 

By  starting  with  the  last  stage  first,  equation  (7)  becomes: 

(    T  N-l  T        1 

J(N)=  minimum  {  xi(N)Qx(N)+  I      x  (k)Qx(k)f 

L—     —     k=i  —     —   J 


or 


[xT(N)Qx(N)j 


J (N) =  minimum Jx  (N) Ox (N)f  +J(N-1)  (14) 

T 

It  will  now  be  noted  that  only  the  x  (N)Qx(N)  term  contains  the  control 
policy  u(N-l);  therefore,  to  minimize  this  function  over  the  last  stage, 
we  need  only  to  find: 


3(xT(N)Qx(N)) 


9  u(N-l) 


=  0  (15) 


T 
By  substitution  of  equation  (4)  into  x  (N)Qx(N) : 


and 


xT(N)Qx(N)  =  r<f>x(N-l)+Au(N-l)lTQ  [*x(N-l)+Au(N-l)J         (16) 


a3Xl(N-U(N)>  •   2aV(N-1)+  2ATQA  u(N-l)=0         (17) 


or 


T 

A  Qj 
•'•  U(N"1>=  "  ~^2  ^(N_1)  (18) 

Equation  (18)  now  represents  the  optimum  control  policy  u(N-l).   It 
will  be  noted  that  we  are  dividing  by  a  multiplication  of  three  matrices; 
this  'is  a  valid  operation  because  a  row  matrix  times  a  square  matrix 
times  a  column  matrix  will  always  yield  a  scalar  number  if  the  multi- 
plication is  valid.   The  multiplication  is  valid  if  the  dimensions  of 
the  respective  matrices  are  proper;  and  in  this  case,  it  will  always  be 
valid. 

By  proceeding  back  one  more  stage,  the  control  policy  u(N-2)  may 
be  optimized. 

J(N)=  minimum  fxT(N)Qx(N)  +xT(N-l)Qx(N- 1) J  +J(N-2)      (19) 
Again  by  direct  substitution,  equation  (19)  becomes: 

J(N)=  minimum}   {(fr-^ft  Q$  }{  4>x(N-2)+Au(N-2)  >  TQ 

C  L   4  Qa.  -I 

[{^^^}{*x(N-2)+Au(N-2)}l 

Ux(N-2)+Au(N-2)    Q  Lx(N-2)+Au(N-2)  I  \ 


(20) 
■Ux(N-2)+Au(N-2)  I  ~Q  |cf>x(N-2)+Au(N-2)  I  \   +J(N-2) 


Minimizing  with  respect  to  u(N-2),  remembering  that  J(N-2)  is  not  a 

=  0: 


function  of  u(N-2);  and  QJ(N-2)/  ,„   nx 

7"£u(N-2) 


^N-2)  =  2ATQ())x(N-2)+2ATQAu(N-2) 


+2AT{^  -  AATQ^}T  *TxT(N-2)Q{*  -*£&■   }A 
ATQA  A  Ok       ~ 


AA^g£    T  AATo<f> 


+2U"i^T}-  u"7^~-u(N_2)=0 


(21) 


Let 


1    L   AT(Q+Pq)A 


<J>,  where  P  =  0 
o 


therefore 


^J(N) 
3.u(N-2)  ■  AxQ(j)x(N-2)+A>^Q  ^Ax(N-2) 


-  J 


T,T, 


(22) 


or 


+ATQAu(N-2)+AT^Q^1A.u(N-2)=  0 

:.  U(N_2)=  _  ±1^1 ayj 


x(N-2) 


(23) 

(24) 


Let 


P  =  ^ 

1    1   1 


(25) 


A  (Q+P  )cf> 

u(N-2)=  -  _,^ 1 —  x(N-2) 

AT(Q+P1)A  " 


(26) 


The  optimum  control  policy  for  the  (N-2)  stage  has  been  found  represent- 
ed by  equation  (26). 

To  further  the  establishment  of  recursion  formulas,  the  next  stage 


back  must  be  considered  and  minimized  with  respect  to  u(N-3).   The  fol- 
lowing symbology  will  be  utilized  in  this  stage  of  the  development: 

A  =  n>x(N-3)+Au(N-3)J  (27) 

and 

B  =  1W4-   L_V|  (28) 

L1       L  \  at(q+p1)a/J 

The  cost  function  then  becomes: 

J(N)  =  minimum  {ATQA+BTQB}+J (N-3)  (29) 


Let  *2  =  [i-  ^l)L 


AT(Q+P1)A 


Therefore 


(30) 


J(N)=  minimum{ATQA+   ti>   i\> A~|    TQ  [^^l     }+J(N-3)  (31) 


Minimizing   the   cost   function  with  respect   to  u(N-3) 


3^3)   =   2ATQ<J,x(N-3)+2ATQAu(N-3) 


+2A%0  TJo>   4>   <J>x(N-3) 
-      2      1      i    2 


+2AT^2%1TQ^1^2-u(N_3)    =   °  (32) 

Solving   for  u(N-3): 

u(N-3)-  -     L-L-U        x(N-3)  (33) 

atq  a+at  ^  i  <W  i  i>  2  A 


' 


giving 


u(N-3) 


AT(Q+AV^J* 
2  ] 1  2 

AT(Q+^Tq^2)A 


x(N-3) 


(34) 


Let 


Therefore 


p2=  *2*iQ*1*2 


U(N_3)  =  -  A  (Q+p2) 


Ai(Q+P2)  A 


x(N-3) 


(35) 


(36) 


It  is  believed  that  an  iterative  relationship  may  now  be  shown  without 
regressing  to  another  stage.   The  following  equations  may  readily  be 
shown  to  be  true  by  an  inductive  proof. 


where 


or 


and 


Jk+1 


m   [x  _  MT(Q+Pk)  7 

AT(Q+Pk)  A* 


Pk   =  A  T\)       ,   P   -   0 

K   rk  k-1  k   ° 

u(N-k)  ,  -  AT(Q^k-i)  *   i(H0 

AT(Q+Pk_1)  A 


(37) 


(38) 


(39) 


(40) 


x(k+l)  =  <J>x(k)+Au(k) 


(41) 


To  illustrate  the  use  of  the  above  equations,  the  system  presented 
in  Chapter  I  will  be  considered  with  the  sampling  period  equal  to  one 


second: 


<J>  = 


1    1 
0   1 


A  = 


Q  = 


l   o 

0   0 
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To  find  u(k),  equations  37,  39,  and  40  must  be  solved: 


and 


and 


*  -  ♦  -  m'q*   = 


ATQA 


0   0 
-2  -1 


u(k)=  - 


Pj_  -  ^Q*x  =  0 


ATQ(j) 

=^s — x(k)  =  -2x1(k)-2x9(k) 

A  QA 


(42) 


(43) 


(44) 


By  substitution  into  equation  (44)  and  solving  for  u(0),  we  find 
x1(l)=0  and  x2(l)=  -2.   The  following  phase  plane  (Figure  II-l)  results 


1- 


t 


~i ? 

2   x. 


.1- 


4 


I  / 

/ 


-2T 


Figure  II-l:   Phase  Plane  For  Example  Problem 

This  clearly  indicates  that  a  limit  cycle  results  and  that  the  cost  func- 
tion, as  determined  with  the  use  of  this  weighting  matrix,  yields  a  rela- 
tive minimum  cost  of  one  unit. 

Using  the  derived  recursion  formulas  37,  39,  40  and  41,  PROGRAM 
0PC0N1  was  written.   Figures  II-2  through  11-10  represent  abbreviated 
flow  charts  of  the  program  and  the  subroutines  used  in  the  program.   It 
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must  be  remembered  that  the  arithmetic  operations  represented  in  the 
flow  charts  are  all  matrix  operations.   As  can  be  seen  from  Figure  II-2, 
the  required  inputs  to  the  program  must  be  read  from  data  cards.   The 
following  table  (Table  II-l)  indicates  the  order  that  the  data  cards  are 
read,  the  format  required  on  these  cards,  and  the  definition  of  the  vari- 
able concerned.   The  primary  output  of  0PC0N1  is  the  value  of  each  state 
variable,  the  control  signal,  and  the  cost  of  each  stage.   A  complete  pro- 
gram listing  may  be  found  in  the  appendix. 


Variable 

Definition 

Format  of 
Data  Card 

M 

Number  of  stages  plus 
one 

15 

N 

Order  of  system 

15 

A 

Column  vector  of 
system  impulse 
responses  due  to 
control  inputs 

8F10.6 
(read  from  top  of 
column  to  bottom) 

<j> 

State  transition 
matrix 

8F10.6 
(read  in  order  by 
columns) 

Q 

Weighting  matrix  of 
the  cost  function 

8F10.6 
(read  in  order  by 
columns) 

X  (0) 

Initial  conditions 
on  the  state  vector 

8F10.6 

(read  from  top  of 

column  to  bottom) 

f 

i 

Table  II-l:   Definition  and  Format  of  Data  Cards  for  Program 
0PC0N1. 


1  2 


START 


Y 

■i 

READ: 
M,N,A,<|>, 
Q,x(0) 

" 

PRINT: 
DATA  CARDS 

oJ+1 


Z(1,I,J)=  Q(I,J) 
ZZ(I,J)=  Z(1,I,J) 


C    j    «    ) 


(    I  .»    ) 


al 


Figure  II-2-'  Flow  Chart  Of  Program  0PC0N1 
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al 

V 


FIND: 
PSI,P,and  (Q+P) 


P(I,J,K)=  PJ(J,K) 

ijj(I,  J,K)=  RR(J,K) 

Z(I+1,J,K)=  ZZ(J,K) 


K  :  N 


Q 


>, 


I  ;  M-l 


> 


a2 


FIGURE  II-2:   (continued) 


1* 


a2 


PRINT : 
P,4>,and(Q+P) 


SM(I,J,k)=  Z(M-I,J,k) 


K  :  N 


J  :  N 


> 


I  :  M 


> 


PRINT: 
INVERTED  (Q+P) 


Figure  II-2  :  (continued) 
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SMM(J,K)=  SM(I,J,K) 


(   K  :  N   J 


GHZ) 


FIND: 


-H 


u(I)=  H*x(I) 


x(I+l)=  <J)*x(I)+A*u(I) 


FIND: 
COST(I) 


C    i  ••  M-Q 


>     K. 


FIND: 
COST(M) 


PRINT: 
u(I),H(I),x(I),COST(I) 


A 

STOP 


Figure  II-2:  (concluded) 
16 


START 


J-C 


i t 


A(I,J)=0.0 


A(I,J)=  A(I,J)+B(I,K)*C(K,J) 


(  K  :  N   ) 


(   J  :  M  ) 


(  i;i  ) 


EXIT 


Figure  II-3:  Flow  Chart  For  Subroutine  PROD(A,B,C,L,M,N) 
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START 


:j+i 


A(I,J)  =  B(I,J)+C(I,J) 


c 


ZD 


C 


o 


A 

EXIT 


Figure  II-4:  Flow  Chart  For  Subroutine  SUM(A,B,C,N) 
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B(J,I)  =  A(I,J) 


c 


J  :  N 


<   ^ 


c 


D 


I  :  N 


3 


EXIT 


ieure  II-5  :  Flow  Chart  For  Subroutine  TRANSQ(A,B,N) 


Figur 
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START 


1 

i 

1       i 

E 

B(1,I)=  A(I,1) 

A 

j..  ,  ,-  .-   ..  .1 

1 

V 

(     I.H       )                          If 

7 

< 

\ 

EXIT 


Figure  II-6:  Flow  Chart  For  Subroutine  TRANCOL(A,B,N) 


START 


FIND: 


±. 


MO  =  Q  *  x 


MON  =  xT*MO 


EXIT 


Figure  II-7:  Flow  Chart  For  Subroutine  Money (MON, E,Q,N) 
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START 


FIND: 

AT 

< 

' 

AA  =  AT*(Q+P) 

»' 

BB  =  A  *  AA 

\' 

CC  =  (Q+P)*  A 

' 

• 

EE  =  AT*CC 

i 

• 

A  =  I 

• 

' 

BB  =  BB/EE 

' 

BB  =  A  -  BB 

• 

\\)   =  BB*<j> 

EXIT 


Figure  II-8:  Flow  Chart  For  Subroutine  PPSI(B,C,D,E,F,N) 
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START 


FIND: 


D   =   Q   *   ty 


P   =   ip      *   D 


EXIT 


Figure  II-9  :  Flow  Chart  For  Subroutine  PP(A,B,C,N) 


START 


E  =    (Q+P)*    0 

1 

' 

FIND: 

AT 

1 

' 

H  =   AT    *  E 

-     —           -"  — ' 

G   =  (Q+P)*  A 


A  =  AT  *  G 


H  =  H/A 


EXIT 


Figure  11-10:  Flow  Chart  For  Subroutine  FF(A,B,C,D,N) 
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CHAPTER  III 
PROGRAM  OPCON2 

Program  0PC0N2  was  written  as  the  final  phase  for  a  versatile  simu- 
lator to  determine  a  constant  gain  matrix  for  optimization  of  a  control 
system  according  to  the  following  cost  function: 

J(N)=minimum  ,  8   (xT(k)Qx(k)+ru2 (k-1 ) )  (45) 

k-i  - 

Again  u(k)  is  a  scalar  control  and  the  Q  matrix  is  determined  in 
the  same  manner  as  outlined  in  the  previous  chapter.   The  variable,  r, 
is  a  weighting  constant , (Lagrange  multiplier)  determined  by  user  needs 
and  acts  as  a  quadratic  integral  constraint  on  the  control  effort;  such 
as  a  limit  on  energy  necessary  to  provide  the  desired  control.   By  let- 
ting r  equal  zero,  we  have  the  cost  function  of  0PC0N1  and  identical 
results  are  obtained  from  the  two  programs.   With  the  use  of  0PC0N1 ,  no 
restraints  are  put  on  the  control  required  by  the  system  which  may  re- 
sult in  the  necessity  of  an  unlimited  supply  of  fuel  to  obtain  the  de- 
sired optimal  results. 

Again  let  us  consider  the  system: 

x  =  Fx  +  Du  (1') 

whose  discrete  solution  is: 

x(k+l)=  <j»x(k)+A_u(k)  (4') 

where  u(k)  and  A_  are  defined  in  the  same  manner  as  in  the  previous  chapter, 

By  starting  with  the  last  stage  first,  equation  45  becomes: 

J(N)=minimumj  xT(N)Qx(N)+ru2(N-l)|  +J(N-1)  (46) 

or  by  substitution  of  equation  4': 

J(N)=minimumf(<j>x(N-l)+Au(N-l))  Q  (<j>x(N-l)+Au(N-l)  )+ru2  (N-l)]  +J(N-1) 

[  J  (47) 
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We  now  minimize  J(N)  with  respect  to  u(N-l)  noting  that  J(N-l)  is  not  a 

function  of  u(N-l). 

3  J(N)         T         T 

3u(N-l)  =  °  =  A^^xCN-D+^QAuCN-D+ruCN-l)  (48) 


or, 


J, 


U(N_1)=  _  J_Q1_   N_  (49) 

AjQA+r  _ 

As  in  Chapter  II,  the  denominator  of  equation  49  will  always  be  a  scalar 
number;  therefore,  no  matrix  inversion  is  required. 

By  proceeding  back  one  more  stage,  the  control  policy  u(N-2)  may  be 
optimized.   Equation  45  becomes: 

J(N)=minimum  xT(N)Qx(N)+xT(N-l)Qx(N-l)+ru2(N- 1)  +ru2 (N-2) |+J(N-2)    (50) 
The  following  symbology  will  be  used  for  clarity  in  this  minimization: 

A  4  <|>x(N-2)+Au(N-2)  (51) 

and, 

B4_^!_  (52, 

A_  QA+r 

By  direct  substitution  using  4',  51,  and  52,  equation  50  becomes: 

J(N)=minimum{(c})A+ABA)TQ((()A+ABA)+ATQA+r(BA)T(BA)+ru2(N-2)j+J(N-2)     (53) 
By  rearranging  terms, 
J(N)=minimum)ATQA+(C(|)+AB)A)TQ(C(})+AB)A)+r(BA)2+ru2(N-2)[+J(N-2)      (54) 


We  now   let, 


*[   =  B  =  -     -* (55) 

A_TQA_+r 


and 


*1-  4>+ha[  (56) 

Equation  54  now  becomes: 

J(N)=minimumjATQA+A%  Qi\i  A+r (a^A) 2+ru    (N-2)] +J (N-2)  (57) 

2  4 


We  now  minimize  J(N)  with  respect  to  u(N-2), 

9   J  (N)  t  t    T  T 

3u(N_2)      =0=   AiQA+Ai4;1QlJ;iA+ra1Aa1A+ru(N-2)  (58) 

Solving  for  u(N-2)  by  substitution  of  51  into  equation  58, 

u(N-2)=  -   /ATQ+ATiJj?Q^1+rATa1a1  f  <J> 

K        '  ~  —       i,       \       .     x(N-2)  (59) 

IT  TT  T  T        T     —  v-»'/ 

JAiQA+Ai4;i1Qi|;1A+rA  a]_a1    A+rj 


We  now   let, 


P1=4;iQ^l+Q+ra1a1  (60) 


Equation  59  becomes: 

u(N-2)=  -  ATM 

ATP1A+r 


x(N-2)  (61) 


As  we  did  in  the  previous  iteration,  we  let, 


T=  _  ATP!* 


and, 


a2     -^ (62) 

ATP  A+r 


4>2=   ())+Aa2T  (63) 


By  proceeding  back  one  more  stage,  a  recursive  relationship  may  be  defined 
We  shall  define  the  following  quantity  for  brevity  in  this  development: 

G£<}>x(N-3)+Au(N-3)  (64) 

The  cost  function,  45,  may  be  written  as, 

J(N)=minimum  xT(N)Qx(N)4-xT(N-l)Qx(N-l)+x  (N-2)Qx(N-2) 

+ru2 (N-l)+ru2 (N-2)+ru2 (N-3)  +J (N-3)  (65) 

By  substitution  of  4',  55,  56,  62,  63,  and  64,  equation  65  becomes: 
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J(N)=minimum  j (^2<hG)TQ(>MlG)+r  (^2aiTG)T(>halTG) 

+  02G)TG(^2G)+r(a2  G)T(a2TG) 

+GTQG+ru2(N-3)j+J(N-3)  (66) 

Minimizing  the  cost  over  N  stages  with  respect  to  u(N-3), 

3  J(N)     _    AT,T,T  ,  .    AT,T  ,    .T 
3u(N-3)  =  °  =  -  Vl^lV*-  M^2G+-  QG 

+rAT^Ta  a^  G+rATa  a^G+ru(N-3)  (67) 

Substituting  64  into  equation  67  and  solving  for  u(N-3),  we  find  the 
optimum  control  policy: 

u(N-3)    = 

{  AT^2^lQ^l^2+-T^2Q^2+-TQ+r-T^2ala1^2+r-Ta2a2)  * 


( kT\bT\bTQ\b   \b   A+ATi|;TQtJj   A+ATQA+rAT^Ta  aT\b   A+rATa   aTA+r]      — 
\-     2    1      1    2 r2      2 ~    ~     2   1    12-    -     2   2-     J 


x(N-3)      (68) 


Combining  terms, 

u(N-3)=  _  {A  ^2(*lQ*l+Q+raiai)^2+A  Q+rA  a2a2  }  * 

I  ATiJjJ(^fQi(;1+Q+ra1a^)ij;2A_+A_  QA+rA  a2a2b+r*- 


x(N-3)      (69) 


From  equation  60,  equation  69  becomes 


„fN-^-   -  f  AT^2Pl^2+ATQ+rATa   a0T}   <{> 

U(N   3)~       '         ^    i    ^  -  h     -     ?~?    i x(N-3)  (70) 

{  AT^2pi^2A+ATQA+rATa2a2TA+r} 

We  now   let, 


P2=^P1^2+Q+ra2a2T  (71) 


and  equation  70  becomes: 


T 

u(N-3)  =  -  -  gaJ   x(N-3)  (72) 

A  P2A+r 


Let, 
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a.T=  -  ATP0* 


ATP  A+r 


and, 


i\i  =   <()+Aa3 


(73) 


(74) 


As  in  Chapter  II,  a  recursive  relationship  may  now  be  shown  without 
proceeding  back  another  stage  and  may  be  shown  to  be  true  by  an  induc- 
tive proof.   Let 

J  =  _  Aipk-i* 

9  T 


and, 


where, 


T  =  _  Alpk-i4> 
k     ATpk-iA+r 


ax  =  0 
o 


\b     =  cb+Aa   .  \b   -  0 
k    -  k  '   o 


P  =  ^TP   \\>   +Q+ra  a  T,  P  =0 
k   k  k-1  k  y   k  k    ° 


(75) 


(76) 


(77) 


and, 


u(N-k)=  a,  x(N-k) 
k^ 


(78) 


from, 

x(k+l)=  <t>x(k)+Au(k)  (79) 

By  letting  r  equal  zero  and  solving  the  same  illustrative  example  as 
in  the  previous  chapter,  the  same  results  are  obtained  as  illustrated 
below. 


and, 


and, 


A   QA+r 


-2.0      -    2.0 


\b   =   <t>+Aan  ' 
1  _  1 


0  0 

-2.0  -1.0 


rp  T* 

Px-  ^1Q^1+Q+ra1a1   =Q= 


1.0        0 
0  0 


(80) 


(81) 


(82) 
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and  therefore, 


u(k)=aiTx(k)=  -  2x  (k)-2x2(k)  (83) 


Equation  83  indicates  the  same  results  as  in  Chapter  II  and  perscribes 
the  same  phase  plane  (Figure  II-l). 

Through  the  use  of  the  recursion  formulas  75,  76,  77,  78,  and  79, 
PROGRAM  0PC0N2  was  written.   Figures  III-l  through  III-5  and  Figures  II-3 
through  II-6  represent  abbreviated  flow  charts  of  the  program  and  the  sub" 
routines  used  in  the  program. 
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START 


V 


READ  AND   PRINT 
M,N,IT,R,A, 
cj)  ,Q,and  X 


KN=  0 
TCOST=  0.0 


t i 


P1(I,J)=  Q(I,J) 


c 


J  :  N 


3 


CUD 


> 


A 

v 


al 


Figure  III-l:   Flow  Chart  for  PROGRAM  OPCON2, 
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al 


FIND: 


AT 


ALFAT(M-I,J)=  AT(1,J) 


J:N 


> 


FIND: 


y\>   and  P 


KN=  KN+1 


KN:IT 


>, 


(PRINT: 

AT,tJ;,  and  P 


KN=  0 


^T 


a2 


h\ 


Figure  III-l:   (continued) 
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P1(J,K)=P(J,K) 


K:N 


J:N 


>, 


I:M-I 


vl- 


AT(1,I)=  ALFAT(1,I) 


I:N 


a3 


61 


Figure   III-l:       (continued) 
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CJ+2 


/ 


AT(J,K)=   0.0 


K:N 


» 


J:N 


>, 


FIND: 


u(l,l) 


Z=   u(l,l) 
L   =    1 


FIND: 

DOL 


TCOST   =   TCOST+DOL 


a4 


Figure  III-l:   (continued) 
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a4 


y 

PRINT: 

L,Z,X(1) ,    and   TCOST 

>r 

KN=    1 

Y(J,1)=  X(J,1) 


c 


J:N 


X=<j)*Y 


DELTA(J,1)=Z*DEL(J,1) 
X(J,1)=  X(J,1)+DELTA(J,1) 
AT(1,J)=ALFAT(I,J) 


J:N 


T 

a5 
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Figure   III-l:       (continued) 


3  3 


AT(J,K)=  0.0 


K:N 


J:N 


FIND: 


u(I,l) 


Z=  u(I,l) 


FIND; 


DOL 


TCOST=  TCOST+DOL 
KN=  KN+1 


a6 
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A 
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Figure  III-l:   (continued) 
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a6 

L 


KN:IT 


32 


PRINT: 

I,Z,X,  and  TCOST 


KN=  0 


I:M 


* 


STOP 


Figure  III-l:   (continued) 
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START 


FIND:  T 

A 


ab=  ax*p 

AC=  AB*A 
AD=  AB*<j) 


AT(1,I)=  -AD(1,I)/(AC(1,1)+R) 


I:N 


A 

EXIT 


Figure  III-2:   Flow  Chart  for  SUBROUTINE  ATRAN(AT,P,PHI,DEL,R,N) 
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START 

V 


AB=   A*AT 
ip=   <J>+AB 


EXIT 


Figure  III-3:   Flow  Chart  for  SUBROUTINE  PPSI(PSI,PHI,DEL,AT,N) 


START 


FIND: 
YT 


A  =  YT*Q 

AA  =  A*Y 

DQL=  AA(1,1)+R*Z*Z 


A 


EXIT 


Figure  III-4:   Flow  Chart  for  SUBROUTINE  COST(DOL, Y,Q,R,Z,N) 
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J+ 


START 


y   I 


AA(I,1)=  AT(1,I) 


I:N 


>/ 


FIND: 


AB=    <^T*P1 
AC=  AB*i> 
AD=  AA*AT 


rj-*-i> 


!! !! 


AD(I,J)=  R*AD(I,J) 


J:N 


>• 


o7 


A 
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Figure  III-5:   Flow  Chart  for  SUBROUTINE  PP(P,PSI,P1 ,Q,AT,R,N) 
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a7 

I 


I:N 


>/ 


I 


< 


AE=  AC+Q 
P=  AE+AD 


A 

EXIT 


63 


Figure   III-5:       (concluded) 


39 


Figure  III-l  shows  that  the  inputs  to  the  program  must  be  read  from 
data  cards.   Table  III-l  shows  the  order  in  which  these  cards  are  read, 
the  format  required  on  each  card,  and  the  definitions  of  the  variables 
concerned. 


Variable 

Definition 

Format  of 
Data  Card 

M 

Number  of  stages 
Plus  one 

19 

N 

Order  of  system 

19 

IT 

Frequency  of  the 
desired  printing 
of  output 

19 

r 

Weighting  constant 
of  cost  function 

F10.6 

A 

Column  vector  of 
system  impulse 
responses  due  to 
control  inputs 

8F10.6 
(read  from  top  to 
bottom  of  column) 

4> 

State  transition 
matrix 

nF10.6 
(read  in  a  row  at 
a  time) 

Q 

Weighting  matrix  of 
the  cost  function 

nF10.6 
(read  in  a  row  at 
a  time) 

x(0) 

Initial  conditions 
on  the  state  vector 

8F10.6 
(read  from  top  to 
bottom  of  column) 

Table  III-l:   Definitions  and  Format  of  Data  Cards  for 

PROGRAM  0PC0N2. 
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The  output  of  PROGRAM  0PC0N2  is  the  value  of  each  state  variable, 
the  control  signal,  and  the  total  cost  of  all  previous  stages  includ- 
ing the  present  stage.   The  amount  of  printed  output  may  be  varied  by 
the  use  of  the  IT  variable's  input  data  card.   If  all  output  is  to  be 
printed,  then  IT  must  be  set  to  one  by  the  data  card.   If  the  user 
wishes  to  have  1,000  stages  and  has  a  requirement  of  print  out  of 
only  a  sampling  of  these  stages,  IT  should  be  set  accordingly.   By  let- 
ting IT  equal  ten,  only  every  tenth  stage  will  be  printed  out.   Another 

T 
valuable  output  is  the  value  of  a  which  approaches  the  constant  gain 

K 

matrix  necessary  for  optimal  control  of  a  continuous  system.   This  is 
achieved  by  using  a  small  sample  period  and  many  iterations  thereby 
approaching  the  simulation  of  a  continuous  system.   0PC0N2  is  limited 
for  this  use  because  the  computer  memory  size  limits  the  program  to  ap- 
proximately 2,000  stages.   If  the  user  has  no  need  for  the  outputs  of 
the  values  of  the  control  signals,  state  variables  and  total  costs, 

then  the  program  is  unlimited  in  its  number  of  stages  because  the  value 

T 
of  a  may  be  printed  out  as  often  as  desired  without  storing  it  in  mem- 
K 

ory.   This  alteration  has  been  made  in  PROGRAM  0PC0N3 ,  and  a  listing 
of  this  program  may  be  found  in  the  appendix. 

The  versatility  of  this  program  may  be  further  shown  in  using  it 
to  simulate  an  optimum  control  in  accordance  with  the  cost  function, 


i["xT(N)X(N)J 


J(N)=  minimum  X1 (N)X(N)  (84) 

which  we  shall  call  the  final  value  criterion.   Using  the  same  method 
of  derivation  as  before,  we  find  the  following  recursion  formulas: 
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ty  -♦+Aai>   4^  =  0 
K        —  K        o 


(85) 


and, 


T  ATP        d>  t      „ 

T 


(86) 


where, 


PTr     =  \b       P        *        ,P  =   I 
K-l      \-l  K-21-l'    ° 


(87) 


and, 


u(N-K)=  aAX(N-K) 
K 


(88) 


To  use  0PC0N2  for  this  cost  function,  we  need  only  to  set  P  equal 
to  the  identity  matrix,  Q  to  the  null  matrix,  and  r  equal  to  zero. 

As  an  example  of  this  cost  function,  we  shall  again  use  the  same 
illustrative  system.   We  find  that: 


T     A- 
ai  =  - 


ATA 


[.,  .,.,] 


(89) 


and, 


and, 


and, 


\\>   »^+Aa1  = 


.8     .4 
-.4    -.2 


?1   =   ^  h= 


T 
a  = 

2 


ATpiA 


.2 


[-1.0   -1.5] 


(90) 


(91) 


(92) 


therefore, 


u(0)  =  -  X1(0)  -  1.5X2(0) 


(93) 
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Given  initial  conditions  of  X  (0)  =  1.0  and  X„(0)  =  0.0,  equation 
93  results  in, 

u(0)  =  -1.0  (94) 

therefore, 

X1(l)  =  X1(0)+0.5u(0)  =  0.5  (95) 

and, 


From  equation  89 , 


X2(l)  =  X2(0)+u(0)  =  -  1.0 


u(l)  -  -.4X1(1)  -  1.2X2(1)  =  1.0 


and, 


and, 


X1(2)  =  X1(l)+X2(l)  +  .5u(l)  =  0.0 


X2(2)  =  X2(l)+u(l)  =  0.0 


(96) 


(97) 


(98) 


(99) 


Figure  III-6  is  the  phase  plane  resulting  from  the  calculated 
controls  and  shows  that  the  final  states  at  N  equal  two  are  definitely 
at  a  minimum  and  the  cost  function  equals  zero. 

*X, 


Figure  III-6:   Phase  Plane  for  Example  Problem, 
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With  the  use  of  0PC0N2  and  0PC0N3,  many  control  systems  may  be 
simulated.   Using  these  two  programs,  a  system  may  be  designed  for  a 
number  of  optimization  criteria  while  greatly  reducing  design  time. 
In  the  next  chapter,  a  typical  system  will  be  considered  assuming 
an  optimization  criterion  and  the  resulting  design  will  be  obtained. 
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CHAPTER  IV 
EXAMPLE  PROBLEM 

To  show  the  use  of  the  programs  that  were  written  and  developed 
in  the  preceding  chapters,  the  following  example  problem  was  chosen. 
The  transfer  function  is  typical  of  a  servo  motor,  the  roll  attitude 
control  of  a  missile,  and  many  others.   Figure  IV- 1  shows  the  open 
loop  plant  to  be  controlled. 


u 


Figure  IV- 1:   Block  Diagram  of  the  Open  Loop  Control  System. 
The  following  state  equation  defines  the  system: 


x  = 


0 

1 

x  + 

0 

0 

-1 

1 

•"■* 

.  . 

(100) 


With  the  use  of  PROGRAM  PHIDEL  and  a  sample  period  of  one  second,  <f> 
and  A  are: 


x(k+l)  = 


1.00   .632 
0.0    .368 


x(k)  + 


.368 
.632 


"(k) 


(101) 


To  show  the  flexibility  of  PROGRAM  0PC0N2,  four  different  cost 
functions  will  be  considered  and  compared.  The  first  cost  function 
will  be  the  terminal  value  function  for  time  optimal  control.   (Case  I) 

J(N)  =  minimum  xT(N)x(N)  (84') 

This  will  then  be  compared  to  the  terminal  value  control  with 
an  energy  control  shown  by  equation  (102).   (Case  II) 
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T        ^  2 
J(N)  =  minimum  x  (N)x(N)+>  u*(K-l)  (102) 

k=l 

The  next  case  of  interest  minimizes  the  sum  of  the  squares  of  the 

states  with  an  energy  limitation.   This  is  represented  by  equation  (45'). 

(Case  III) 

k=N 
J(N)  =  minimum  ^  xT(K)Qx(K)+ru2 (K-l)  (45') 

k=l 
where, 

Q  =  I  and   r  =  1.0  (103) 

The  final  cost  function  will  also  be  that  of  minimizing  the  states 
with  an  energy  limitation;  however,  by  use  of  PROGRAM  0PC0N3,  a  steady 
state  feedback  matrix  was  found.   These  values  were  then  used  in 
PROGRAM  0PC0N2  to  obtain  a  transient  solution.   (Case  IV)   The  elements 
of  this  control  are: 

aT(l,l)  =  -0.627797  (104) 

and, 

aT(l,2)  =  -0.615037  (105) 

With  the  use  of  PROGRAM  0PC0N2,  a  transient  solution  was  obtained 
in  each  of  the  four  cases.   Case  II  was  computed  for  a  twenty  second 
trajectory  in  the  phase  plane,  Figure  IV-3.   Cases  III  and  IV  were  each 
run  for  ten  second  trajectories  while  Case  I  was  optimized  in  two  samples, 

Figure  IV-2  is  a  plot  of  control  versus  time  for  each  of  the  cases. 
It  is  of  interest  to  note  that  both  cases  III  and  IV  utilize  the  same 
controls  during  the  first  seventy  percent  of  their  trajectories;  only 
when  the  control  is  very  small  in  the  last  couple  of  sampling  points 
does  the  control  required  differ.   This  thereby  leads  to  identical  trajec- 
tories in  the  phase  plane  indicating  that  constant  gain  amplifiers  in  the 
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feedback  loop  are  just  as  effective  as  variable  gain  amplifiers  for  com- 
pensation of  this  plant.   Case  II  requires  nearly  a  constant  control  ef- 
fort after  the  initial  control  jump  has  been  made.   The  required  control 
for  the  twenty  second  trajector  is  approximately  equal  to  one  half  of 
the  necessary  control  for  a  ten  second  trajectory. 

Figure  IV- 3  is  a  phase  plane  plot  of  the  four  cases.   One  thing  to 
note  here  is  the  trajectory  of  case  II.   As  the  total  number  of  sample 
periods  is  increased,  the  necessary  control  is  decreased  as  indicated 
above.   This  results  in  the  trajectory  approaching  the  x-.  axis  as  the 
number  of  sample  points  increases  indefinitely.   This  means  that  the 
cost  function  used  is  optimized  when  the  final  value  of  the  states  is 
zero  and  the  control  is  nearly  zero.   The  point  in  the  phase  plane  would 
drift  in  requiring  only  a  very  slight  amount  of  control  to  stop  it  at  its 
steady  state  value. 
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Figure  IV-2:   Plot  of  Control  Versus  Time  for  the  Four  Cases. 
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